Ground-state phase diagram of the quantum Ji — J2 model on the square lattice 
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We study the ground-state phase diagram of the quantum J\ — J2 model on the square lattice by 
means of an entangled-plaquette variational ansatz. In the range < J2/J1 < 1, we find classical 
magnetic order of Neel and coUinear type, for J2/J1 < 0.5, and J2/J1 > 0.6 respectively. For 
intermediate values of J2/J1 the ground state is a spin liquid (i.e., paramagnetic with no valence 
bond crystalline order). Our estimates of the entanglement entropy show that such a spin liquid is 
topological. 
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I. INTRODUCTION 

The spin-1/2 antiferromagnetic (AF) Heisenberg 
Hamiltonian with additional AF next-nearest neighbor 
coupling is one of the most studied models, known as the 
Ji — J2, in theoretical physics; 



iJ = Ji 



/ . Si • Sj + J2 2^ S^■ Sj. 



(1) 
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Here S,; is a spin-1/2 operator associated to the ith lattice 
site and the first (second) summation runs over nearest 
(next-nearest) neighbor sitesii 

Model ([T]) is of relevance for experimentally accessi- 
ble compounds?^ and constitutes a canonical example 
of spin system featuring frustration induced by compet- 
ing AF interactions. It is known that, on the square lat- 
tice, the ground state (GS) of ([1]) displays AF long-range 
order of Neel [ordering wave vector q = (tt, tt)] type at 
J2 = (i.e., in the unfrustrated case);^ the Neel order re- 
mains for sufficiently small values of J2 while at large J2 
the AF order is collinear [q = (7r,0) or (0,7r)]<^ To gain 
insight into the nature of the GS at intermediate values 
of J2, several numerical techniques have been applied. 
The exact diagonalization (ED) of the Hamiltonian ([1]) is 
possible only for small lattice size (up to A^ ^ 40 spins), 
therefore, although this approach may give fundamen- 
tal suggestions, a reliable extrapolation of the relevant 
physical observables to the thermodynamic limit is basi- 
cally not achievable. On the other hand, quantum Monte 
Carlo (QMG) schemes based on the imaginary time evo- 
lution of a trial wave function (WF) are essentially exact, 
even for very large lattices, at J2 = 0, being affected by 
the sign problem at finite J2 (i.e., in presence of frus- 
tration). Variational Monte Garlo (VMC) calculations 
overcome the major limitations of ED and QMG being 
not restricted to small lattice size and, by definition, free 
of any sign instability. However, their accuracy depends 
on the guess for the GS WF. Hence, understanding and 
characterizing the GS properties of the Ji — J2 model 
in the maximally frustrated regime (i.e., J2 ^ 0.5 for a 
square geometry) is as difficult as it is interesting. 

ED findings are consistent with a paramagnetic valence 
bond crystalline (VBG) phase intervening between the 



Neel and collinear one<^ This conclusion also arises (al- 
beit with possibly different VBG patterns) from recent 
projected-entangled pair states?^ series expansion and 
spin-wave theory^i coupled-cluster^ hierarchical mean 
field)^ and renormalizcd tensor networki^ investigations, 
whose main objective has been that of elucidating the 
type of phase transition occurring at the boundary be- 
tween the Neel and paramagnetic VBG state. 

Recent VMG calculations pointed out the existence of 
a spin liquid (i.e., paramagnetic with no VBG order) GS 
for model (jlj on a honeycomb lattice ;^^'^^ revitalizing 
the long-lasting debate concerning the possibility of ob- 
serving a similar phase on the square latticeJ^^— 

In this paper we investigate the GS phase diagram of 
dl]) using a VMG approach based on the variational fam- 
ily of entangled-plaquette states (EPS) J^ The EPS WF 
is a general ansatz which has been applied to a variety 
of unfrustrated (bosonic) or frustrated (fermionic) lattice 
problems yielding results of accuracy at least comparable 
to (or better than) that obtainable with different numer- 
ical schemes or alternative variational WF'si ^^i^^'^^ In 
particular it has been employed with success to determine 
the GS phase diagram of ^ on the honeycomb lattice,— 
giving predictions quantitatively and qualitatively more 
accurate than those affordable using different ansatzsjii 
The main quality of the EPS WF is its systematic im- 
provability which allows one to compute increasingly ac- 
curate estimates of physical observables by sequentially 
increasing the plaquette size (i.e., the number of spins 
comprised in a single plaquette) /. 

Wc estimate GS energies and relevant order param- 
eters, i.e., square sublattice magnetization (SSM) and 
various VBG structure factors, in the range < J2 < 1. 
For each system size, up to A^ = L x L = 256 spins 
with periodic boundary conditions (PBG), the EPS WF 
is independently optimized, starting with N square pla- 
quettes of four sites. The size of our plaquettes is then 
consistently and systematically increased up to sixteen 
spins by adding 2Z^/^ + l sites each time that I is changed. 
We extrapolate the estimates obtained on finite lattices 
of different sizes, for a given plaquette size, to the ther- 
modinamic limit and, finally, asses the dependence of our 
results on /. 

Our main findings are the following: a) the system is 



Nccl ordered up to J2 ^ 0.45, while for J2 > 0.6 the 
magnetic order is collinear; b) for intermediate values of 
J2 all the order parameters considered here vanish in the 
thermodynamic limit signaling the emergence of a spin 
liquid; c) by estimating the entanglement entropyi^^— 
we characterize such a spin liquid as topological. 



II. THE ENTANGLED-PLAQUETTE 
VARIATIONAL WAVE FUNCTION 

In this section we briefly recall the basic aspects of the 
EPS WF, referring the reader to Refs. [la, and [13 for 
further details. 

A generic WF for a spin-1/2 lattice system, e.g., the 
one described in model ([T]), can be expressed as 



\^)^J2^in) 



(2) 



where |n) = |ni, n2, . . . , hn)-, rn is the eigenvalue of erf 
and W(n) is the weight of a configuration of the system. 
The EPS ansatz is based on the following idea: cover 
the lattice with N plaquettes of I sites, and identify the 
weight of a given global configuration with the product 
of TV variational coefficients in one-to-one correspondence 
with the particular plaquettc configuration. Hence, 
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where the configuration of the Pth plaquettc is given by 
the values np = ni_p,n2,p, . . . ,n;,p at its I sites. The 
above choice of T4^(n) fully specifies the EPS WF. Such an 
ansatz, when overlapping (i.e., entangled) plaquettes are 
considered (i.e., I > 1) naturally incorporates spin-spin 
correlations, a crucial ingredient to describe the physics 
of a quantum GS. As already mentioned, the fundamen- 
tal property of this variational choice is its systematic im- 
provability obtainable by increasing I. Given its simple 
form, the EPS WF is efficiently optimizable on a com- 
puter; moreover, once optimized, it allows for an evalua- 
tion of the physical observables efficient as well. Both the 
WF optimization and the estimates of physical observ- 
ables are achieved via the Monte Carlo method i^^'^^i^^ 



III. RESULTS 

Figure [T] shows the GS energy per site e = E/N 
as a function of the system size for various values of 
J2 and different plaquettc size (left part). Regardless 
the system size and the value of J2 our estimates sen- 
sibly improve when / increases. For example we ob- 
tain e{N = 100,1 = 4, J2 = 0) = -0.65900(6) and 
e{N = 100,1 = 16, J2 = 0) = -0.67000(5) using 
2x2 and 4x4 square plaquettes respectively. Esti- 
mates obtained with the same I have been extrapolated 
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FIG. 1: (color on line). Left: GS energy per site of the Ji — J2 
model on the square lattice with PBC, as a function of the 
system size. Estimates are obtained with A'^ plaquettes of 
I = 4 (boxes), I = 9 (triangles) and I = 16 (circles) sites. 
Each dashed line is a fit to numerical data with the same I 
value. Right: GS energy per site, extrapolated to the ther- 
modynamic limit, as a function of the plaquette size. Each 
dashed line is a fitting function linear in Z"'^" (see text). The 
QMC estimate in the J2 = case is also shown for comparison 
(star).^ 



to the thermodynamic limit assuming the scaling form 
e(7V, Z, J2) = e(oo,Z, J2) + a{l,J2)N-^/'^. It is clear that 
extrapolated energies depend on I. While in principle 
larger plaquette sizes have to be employed in order to 
observe numerical convergence of the results, available 
computational resources limit the largest plaquette size 
which can be used in practice. However, by construc- 
tion of the EPS ansatz, the estimate of any observable 
must approach the exact GS value in the limit of large 
I. Therefore, the point here is whether it is possible to 
give a reliable estimate of such a limit based on the re- 
sults obtained with plaquettes of relatively small size. As 
pointed out in a previous worki^ on model ([T]) on a hon- 
eycomb lattice, this is indeed the case and estimates of 
the energy and all the relevant physical observables can 
be obtained, in the large I limit, by those computed with 
plaquette of different size (up to Z = 16 in this study). 
In the right part of Fig. [1] we present an example of this 
procedure: we expand our energy per site extrapolated 
to the thermodynamic limit in powers of l^^/^ and fit 
our data with the smallest number of parameters. We 
observe that a good description of our results is obtained 
employing a single power (i.e., l~'^^^). For instance, in 
absence of frustration, our energies e(oo, I, J2 = 0) fall 
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FIG. 2: (color on line). Left; Neel SSM of the Ji- J2 modelon 
the square lattice with PBC, as a function of the system size. 
Estimates are obtained with A'^ plaquettes of Z = 4 (boxes), 
I = 9 (triangles) and / = 16 (circles) sites. Each dashed line 
is a fit to numerical data with the same / value (see text). 
Right: Neel SSM, extrapolated to the thermodynamic hmit, 
as a function of the plaquette size. Lines are functions bult 
to infer the I dependence of our data (see text). The QMC 
estimate in the J2 = Q case is also shown for comparison 
(star).— Inset: pictorial representation of classical Neel order 
where "up" ("down") spins reside on sublattice a (/3). 



almost perfectly on a straight line when plotted against 
Z""^/^, and we find an extrapolated to the large I limit 
energy in agreement with the QMC estimate;^ exact in 
this case.— 

Aimed at identifying extended regions of the phase di- 
agram which are magnetically ordered (disordered), we 
compute the SSM defined as 



xn^{N) 



1 

iV2 






je/3 



(4) 



where the two summations run over lattice sites belong- 
ing to different sublatticcs. In this work we consider two 
types of magnetic long-range order, schematically illus- 
trated in the inset of Fig. [5] (Neel type) and Fig. [3] 
(coUinear type). Specifically, if the sublattices a and /3 
are chosen according to the inset of Fig. [5] (Fig. [3|), 
Eq. (d]) defines the Neel (collinear) SSM, or the mag- 
netic structure factor at ordering wave vector q = (tt, tt) 
(q=(7r,0)) 

The Neel SSM as a function of the system size, for 
different I and various J2 values is presented in Fig. [51 
For each plaquette size the extrapolation to the ther- 
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FIG. 3: (color on line). Same of Fig. [2] where the SSM 
describes now magnetic order of collinear type (sublattices a 
and /3 are chosen according to the inset). 



modynac limit has been performed assuming the func- 
tional form m2(7V, Z, J2) = m2(cx), I, J2)+/3{1, J2)N-^/^ + 
S{1, J2)N~^. The right part of the figure shows the ex- 
trapolated values m^ {00 , 1 , J2) versus l~^/^. At J2 = 
we find that the Neel SSM, in the large / limit, extracted 
by fitting the data with a second order (in Z~^'^) polyno- 
mial, is in agreement with the QMC result)^ For larger 
J2 this order parameter decreases. At J2 = 0.45 it is 
still finite; this conclusion is based on the observation 
that both a linear (dotted line) and quadratic (dashed 
line) fit of the m^(oo, I, J2 = 0.45) values are acceptable 
and yield a finite extrapolated value in the large / limit. 
Specifically, when two powers are considered, the visi- 
ble upward bending of the data is better described, and 
the extrapolated value of the order parameter is slightly 
larger than that found with a linear fitting function. At 
J2 = 0.5 a fit which includes only the zero-th and first 
order term (in l^^^^) gives an unphysical extrapolation, 
while the inclusion of an extra term, proportional to Z~^, 
leads to a null (within the accuracy of our calculation) 
value of the Neel SSM in the large plaquette size limit. 
Based on this analysis it is reasonable to locate the phase 
boundary corresponding to vanishing Neel order in the 
vicinity of J2 = 0.5. 

For 0.5 ^ J2 ^ 0.6 no evidence of magnetic order (ei- 
ther of the Neel or collinear type) is found, while for 
larger J2, up to J2 = 1 (i.e., the maximum J2 value stud- 
ied here) the system is collinearly ordered. 

The dependence of the collinear SSM on the system 
size is shown, for J2 ~ 0.65, in Fig. [3] (left part). Our 
data extrapolated to the thermodynamic limit versus 
l~^'^ are presented in the right part of the figure. In the 
limit of large plaquette size the collinear order parameter 
stays clearly finite. Such a conclusion is consistent with 
both a linear and quadratic fit of the data. 

Our prediction of an extended paramagnetic region in 
the phase diagram of the Ji — J2 model on the square 
lattice is in agreement with previous works carried on 
with different numerical approaches.— "—'^^'i^ Most of 
them^"— found that VBC order, either of plaquette- or 
columnar-dimer type, characterizes such a magnetically 
disordered region. 



In order to elucidate the existence of VBC order we 
compute, in the parameter range of interest, both the 
plaquette VBC and columnar-dimer structure factor de- 
fined as in Ref. |a. Results for the plaquette VBC struc- 
ture factor are presented, as a function of the system 
size, for two value of J2 at which magnetic order is not 
present, in Fig. 2] The numerical data for finite sys- 
tem size are extrapolated to the thermodynamic limit by 
means of a fitting function quadratic in the inverse of the 
system size. We find that, in the limit of infinite system 
size, the order parameter of our interest vanishes, within 
the accuracy of our calculation, regardless the plaquette 
size used in the EPS WF (the error on the extrapolated 
estimates is of the order of 10""'). The columnar-dimer 
order parameter vanishes as well in the thermodynamic 
limit, being, for finite N, considerably smaller than the 
plaquette VBC one. 

In a recent investigation, based on a renormalized ten- 
sor network calculation, plaquette-like order has been 
found in the thermodynamic limitJ^ We note that our 
GS energies are lower than those obtained in Ref. [13. For 
example, for a system comprising 256 spins at J2 = 0.5 
we find, with the simplest and economical EPS WF, 
based on 2 x 2 entangled plaquettes, e^^^{N = 256, 1 = 
4, J2 = 0.5) = -0.46299(3) being e^^^{N = 256, J2 = 
0.5) = -0.45062 the value reported in Ref. M- 

A VBC phase has also been obtained in Ref. y. Such 
an ordered phase is a consequence of to the relatively 
small bond dimension of the projected-entangled pair 
states adopted, in the mentioned work, to describe the 
the WF of the system, and is found to vanish when a 
larger bond dimension is usedJ^ 

In Ref. 0, the authors cannot rule out the existence 
of a plaquette VBC phase, due to the difficulty of the 
extrapolation to the thermodynamic limit of the relevant 
order parameter. The study presented in Ref. [^ is based 
on the ED, in the subspace of short-range valence bond 
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FIG. 4: (color on line). Finite size scaling of the plaquette 
VBC order parameter for two values of J2 at which the system 
is found in a magnetically disordered GS. Dashed lines are fits 
to the numerical data with the same J2. This order parameter 
clearly vanishes, for all the I values considered in this work, 
in the thermodynamic limit. Data shown refer to the I — 16 
case. 



states, of the Hamiltonian Q. This technique is applica- 
ble, due to the unfavorable scaling with the systems size 
of the computational cost, only to small lattices (i.e., 
comprising up to A^ '~ 40 spins). It is therefore not sur- 
prising that, in Ref. 161, where results for finite lattices 
of maximum size iV = 40 are provided, an accurate and 
reliable extrapolation to the thermodynamic limit of the 
VBC order parameter is not achievable. Conversely, in 
this work, we use a general variational WF (i.e., the EPS 
one) which allows to study the GS properties of lattices 
whose size is considerable larger than that presently ac- 
cessible to ED approaches. This is a crucial aspect, as it 
renders the extrapolation to the large A^ limit of physical 
observables more accurate and reliable. 

Having found that for 0.5 < J2 ^ 0.6 the GS of model 
([1]) is a spin liquid, namely a paramagnet with no VBC 
order, we investigate the nature of such a spin liquid 
by computing the Renyi entanglement entropy ^2. This 
quantity can be defined in a form suitable to Monte Carlo 
estimations^ 
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FIG. 5: (color on line). Upper-left part: dependence of the 
Monte Carlo average of the swap operator on the plaquette 
size for square subregion A comprising a x a spins. Differ- 
ent symbols refer to different a values (a = 4 to 8 from top 
to bottom). Each dashed line is a quadratic (in l^^''^) fit to 
numerical data with the same a. Upper-right part: extrapo- 
lated (to the I — >■ N limit) values of the entanglement entropy 
as a function of the boundary length x of subregion A (di- 
amonds). Also shown, for comparison, the QMC estimates 
(crosses)i^ Lower part: entanglement entropy for two values 
of J2 at which the GS is a spin liquid. Each dashed line is 
a linear fit (to numerical data with the same J2) performed 
to extract the topological entropy 7 (see text). Here yl is a 
I//2 X L subregion of a L x L lattice with PBC. Hence, A has 
no corners. 



ai ,bi ,a2,b2 

= {SWPa) 



Ty(a2,bi)W^(ai,b2) 
V(ai,bi)W^(a2,b2) 

(5) 



where ^ is a portion of the system and B its com- 
plement, B.i (bi), with i ~ 1,2 labehng two copies of 
the system, denotes the spin configuration (e.g., along 
the z-axis) of the lattice sites belonging to subregion 
A (B), W is the weight of a global configuration, and 

SWPa = :i:j:;:^;i:;;!:::^-i is the swap operator^ 

Figure [S] shows the Monte Carlo average of the swap 
operator, for a system of size N = 256 at J2 = 0, as a 
function of the plaquette sizes (upper- left part). Different 
symbols refer to square (a x a) subregion A comprising 
a different number of lattice sites. The I -^ N limit of 
(SWPa), obtained by fitting data with the same a us- 
ing a second order polynomial in l~-^'^ (dashed lines), 
defines, via Eq. ([5|), our estimates of the entanglement 
entropy (shown in the upper-right part of the same fig- 
ure, as a function of the boundary length x of region A). 
Our results for ^2 are in excellent agreement with ex- 
act QMC ones (crosses) j^^ Using the same procedure we 
compute S2{x) at J2 = 0.5 and J2 = 0.55 (lower part of 
Fig. [5]), where the GS of the system is a spin liquid. In 
this case the subregion A is a L/2 x L portion of a L x L 
system, therefore, given that we are using PBC, it has 
no corners. A fit of our data assuming the scaling func- 
tion 5*2 (x) — dx — J allows us to identify the topological 
entropy 7 which is finite (see figure) in both the cases 
studied here. We note that, at J2 = 0.5, our estimated 
topological entropy, 7 = 1.0(2), is consistent (taking into 
account the statistical uncertainty) with that of a Z2 spin 
liquid)^ and in quantitative agreement with that com- 
puted in recent studies based on a density matrix renor- 
malization group calculation on long cylinders fii and on 
a tensor product approach using a cluster update imagi- 
nary time evolution methodiii^ 



of an entangled-plaquette variational ansatz. Values of 
ground-state energy and relevant order parameters have 
been computed for < J2 < 1. In addition to mag- 
netically ordered ground states occurring at J2 ^ 0.5 
(Neel type), and J2 > 0.6 (coUinear type), we find that 
the ground state, in the intermediate region of the phase 
diagram, is a topological spin liquid. Further investiga- 
tions, which are beyond the purpose of this study, are 
needed to understand the nature of the phase transition 
occurring between the Neel phase and the spin liquid 
one as well as between the latter and the collinear phase. 
For example, the Neel to spin liquid appears to be a 
continuous phase transition. However, its full character- 
ization is highly non trivial and a subject of great cur- 
rent interest.— We want to emphasize that the numerical 
approach employed here allowed us to provide accurate 
quantitative predictions which, plausibly in the near fu- 
ture, might be tested experimentally by using ultracold 
atoms in an optical lattice, essentially as quantum sim- 
ulators of the Ji — J2 Hamiltoniani^ Finally, it has to 
be mentioned that the entangled-plaquettte variational 
wave function (as well the methodology employed in our 
study) can be applied, without substantial modification, 
to basically any lattice model. It gives a straightforward 
access to crucial observables such as magnetic and non 
magnetic order parameters or the Renyi entanglement en- 
tropy, constituting, in our view, one of the most reliable 
variational choices to investigate the ground-state physics 
of frustrated (fermionic) Hamiltonians, intractable, due 
to the sign problem, with exact Quantum Monte Carlo 
techniques. 
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